DFT-based calculation of Coulomb blockade in molecular junction 
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Quantum transport through single molecules is very sensitive to the strength of the molecule- 
electrode contact. When a molecular junction weakly coupled to external electrodes, charging 
effects do play an important role (Coulomb blockade regime). In this regime, the non-equilibrium 
Green function is usually substituted with master equation approaches, which prevents the density 
functional theory from describing Coulomb blockade in non-equilibrium case. Last year, we proposed 
an Ansatz to combine the non-equilibrium Green function technique with the equation of motion 
method. With help of it. Coulomb blockade was obtained by non-equilibrium Green function, and 
completely agrees with the master equation results [Phys. Rev. B 76, 045408 (2007)]. Here, by 
the Ansatz, we show a new way to introduce Coulomb blockade correction to DFT calculation in 
non-equilibrium case. And the characteristics of Coulomb blockade are obtained in the calculation 
of a toy molecule correctly. 

PACS numbers: 73.23.Hk, 73.63.-b, 72.10.-d, 85.65. -|-h 
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I. INTRODUCTION 

Single molecule electronicsii^i^ has been mostly inves- 
tigated in the high temperature and strong contact to 
the electrode regime. The opposite limit of low tem- 
perature and weakly coupled molecular junctions poses 
a challenge to the currently available experimental tech- 
niques. Still the possibility to probe the spectroscopy 
of single molecule junctions via a lateral gate could offer 
new insights to the peculiar coupling of the electrical and 
mechanical degrees of freedom at the nanoscale. In order 
to be able to establish the transport mechanisms govern- 
ing such molecular junctions, a technique which could 
tackle on one hand single electron charging effects and, 
on the other hand, the inclusion of the electron-vibron 
coupling is of extreme importance. 

In the last ten years, the nonequilibrium Green func- 
tion (NEGF) formalism has been successfully employed 
to describe transport observables on the base of a den- 
sity functional theory (DFT) description of the elec- 
tronic structure^i^i^i^iiiSiiSiiiSiiiii^ and model Hamiltonian 
approaches r^iiiii^ Recently it is applied for the influence 
of the vibron dynamics onto a molecular transistor, and 
lots of excellent results are obtained fi^ii^i^ 

However, when coming to the CB regime, the NEGF 
method has been usually substituted with master equa- 
tion approaches (ME), which prevents DFT from describ- 
ing CB effects. Last year, we proposed an Ansatz to com- 
bine NEGF with equation of motion (EDM) methodJ^ 
With help of the Ansatz above, in non-equilibrium case, 
the Coulomb blockade can be completely described just 
within single-particle space simply, and the result fully 
agrees with the one from ME which is performed in many- 
particle space. 

For DFT to describe Coulomb blockade,—!^ the dou- 
ble subspace (spin-up and spin-down) is usually em- 
ployed. However, it is hard to obtain the Coulomb block- 
ade effects correctly in non-equilibrium case. Our pur- 
pose is to introduce Coulomb blockade effect to non- 



equilibrium DFT calculation by this Ansatz. 

In this paper, with the help of the Ansatz in Ref J^, by 
the model Hamiltonian and EOM approach, we propose a 
self energy to describe CB effects in non-equilibrium case. 
Electronic occupation number, electronic current and 
differential conductance are calculated. The charging- 
induced steps and the ratio of 2/3 : 1/3 in the step heights 
of the occupation number and the current, which are the 
important characteristics of CB by MEj^ is obtained cor- 
rectly. The comparison with the complete CB resultsi^ 
also is done. For the occupation number and the current, 
the difference is very small. For the conductance, only in 
the peak height, the difference is very clear, while there 
is no difference in the peak position. The CB stability 
diagram is also shown in the paper. We can see that the 
self energy can describe CB characteristics. The more 
important is that it is very convenient to introduce this 
self energy to DFT code even in non-equilibrium case. A 
scheme to perform CB correction in DFT calculation is 
suggested with double counting correction (DCC). Then 
it is realized in gDETB^i^ii^ A toy molecule is taken 
for testing, and the CB characteristics are shown in the 
results correctly. 

The paper is organized as follows: firstly a self en- 
ergy for CB is proposed by the model Hamiltonian and 
EOM approach (Sec. II); secondly, based on the self en- 
ergy above, a scheme with DCC is proposed to introduce 



LI 



CI 



D 



C2 L2 

















Au 










Au 




c c 


— s;^ 








Gate 



FIG. 1: (Color online) The molecule for Coulomb blockade 
calculation. See text for details. 
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CB correction to DFT code (Sec. IIIA); finally, the cal- 
culation on a toy molecule is performed in weak coupling 
regime, and the CB characteristics are shown in the re- 
sults correctly (Sec. IIIB). 

II. METHOD AND FORMULA 

In non-equilibrium DFT calculation here, the CB cor- 
rection will be introduced to each level of the molecular 
fragment D (sec Fig. 1).— Therefore, the multi-level An- 
derson impurity model is read as follows, 

H = H^ + Y,{Hc. + H^ii), (1) 

a 

with 

ra.a m^a 
k.(T 

Ha,D= X! (^".fc,™.'^^l'T,a'^m,^ + 'i-C-)> (4) 

where d and c are the operators for electrons on the dot 
and on the left (a = L) and the right (a = R) lead, Um 
is the charging energy of level m, e^.o- is the (m, a) level 
of the quantum dot, while e^^ is the spin a level of lead 
a in fc space, a =t, i. With the help of the EOM and the 
truncation approximation, we can obtain a closed set of 
equations for the retarded and advanced GFs Gm°a-.n,T, 

/ \^{oi)r / a\f^r / a r s: 

-\~UmG[^J.j^,^^ (5) 
- em,<T - U„, ± iry)G^|;/°^ = {nrn,a)5m,n5a,r 

+E('^);/'^(n™,,)G;„/^^„,,, (6) 

where 

^rn,(T;n,T (('^ni,(T Mn.r)) ^ i C^) 
Q(2)r/a n , | \\r/a (o\ 

and 

W 1 2 

^m,a \^) — ^m.a + ^m,a ~ a i ■ \^ > 

a.k ^^'^ ' 

are the electron self-energies from leads with 77 ~ 0+. 

Re-arranging Eqs. ^ and we can obtain the re- 
tarded GF as follows, 

with 

Q{U)r ^ Q(0)r Q(0)r^i^Q(l)r ^ (11) 

G(°)'- = {lo-H^+ ir;}-^ , (12) 
G(i)'- = {a'-i/o + f/ + i??}~\ (13) 
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FIG. 2: (Color online) (a) the comparison on the occupation 
number of electrons as a function of bias voltage, (b) the 
comparison on the current as a function of bias voltage. = 
el = -0.5 eV, Vg = 1.0 V, [/ = 1.0 eV, F = 0.05 eV, VL = 
— Va = Vbias/2. (n) = (nj) + The red dash curves are 

from approximation in Eq. [51 the black solid curves from the 
truncation in Ref.— which fully agree with the results from 
ME,— while the blue dot curves are for the case (7 = eV in 
which there is no CB. In the red dash curves and the black 
solid curves, the new steps appear for charging effect, and the 
heights of them are in the ratio of 2/3 : 1/3, which are the 
important characteristics of CB obtained by ME.— 



i/o is a diagonal matrix composed of ^h(= ^h) is 
the one of {nm,s)Um, and U is of Um, while is the 

self-energy matrix from Eq. 

From Eq. PH)) . we can see that the Coulomb interac- 
tion is just included in G^'^-'''. Therefore, with the help 
of Eq. (ini), the retarded CB self energy I](CB)r can be 
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FIG. 3: (Color online) The comparison on conductance with 
the same parameters in Fig. 1. The red dash curve is from 
approximation in Eq. |6] while the black solid curve from the 
truncation in Ref.— which fully agrees with the one from 



obtained by the relation 



[U! — Hq + IT] U! - Hq + IT] 



(14) 
(15) 



uj — Ho — U + i?7 , 

The result is as follows, 

E^CB)'- = - //o - [G^^^l-^ 

By the Ansatz in Refii^, the lesser GF can be written 
out directly, 

with E(")< - iEoTa/aC^), and = i(S^ - S^), 
fai'^) = /(^ ~ while YfJ"" are the diaganol ma- 
trix composed of Sm^</°, and / is the equilibrium Fermi 
function. After the re-arrangement by the way in Refj^, 
we can obtain, 



G< 



G''I]<G", 



(17) 
(18) 

By the helps of I](*^^)< (see appendix [A|). we could 
get S< = E(")< if S(")< ^ 0. 

Therefore, here, the current / can be calculated simply 
by the Landaur formula , ^^'^^ 

I^^JdiuTr {FlG'^FrG'^} • Ihiuj) - Muj)] , (19) 




■1 1 



FIG. 4: (Color online) The CB stability diagram (the contour 
plot of the differential conductance) calculated by approxima- 
tion in Eq.H with = -0.65 eV, = -0.45 eV, U = 1.0 eV, 
r = 0.02 eV. 



instead of the complicated formula in Ref J^i^. 

In the case of double levels (m = 1 and e^ -^ = j^), the 
numerical calculation is performed, which contributes the 
direct comparison between the truncation in Eq. ^ and 
the one in Reff^^ by which CB results fully in agreement 
with the ones from ME can be obtained.— The compar- 
isons of the electronic occupation number (n) and the 
current J as a function of the bias voltage Vbias at fixed 
gate voltage Vg are shown in Fig. 2. Firstly, the dif- 
ference of the results by the two methods is very small. 
Secondly, the CB characteristics are very clear in them: 
the steps appear for charging- induced level-split, and the 
step heights are in the ratio of 2/3 : 1/3, which are 
also obtained by ME as the important CB characteris- 
tics!^ Only in the comparison of differential conductance 
[Q z= dl/dVhias) (see Fig. 3), the difference in the heights 
of the peaks (f/-axis) is very clear, while no difference ap- 
pears in the positions of the peaks (Vbias-axis). From the 
CB stability diagram in the case of m = 1 and | 7^ e? ; 
(Fig. 4), it can be seen that in this case, the approxi- 
mation in Eq. ^ almost includes all CB characteristics 
(the complete CB stability diagram by NEGF is shown 
in Ref-i^). 

It should be noted that although the approximation 
above can just include some CB characteristics cor- 
rectly;^ it is very convenient to be introduced to DFT 
calculation (which will be shown in the next chapter). 



III. SCHEME FOR DFT 
A. scheme 

For CB calculation, the system is partitioned as fol- 
lows (shown in Fig. 1): the molecular fragment D is 
the CB part, the fragments CI and C2 are contacting 
area, and the LI and L2 are the leads. Then we just 
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introduce CB correction in fragment D, while the non- 
equilibrium calculation should be performed within the 
fragment CI - D - C2. 

In DFT, the KS equation H^^* = S'E'* can be re- 
written as, 



(20) 



with S is the overlap matrix. Then an effective Hamilto- 
nian matri x^^'^^ for model-Hamiltonian calculation can 
be obtained from KS one. 



KS 



(21) 



After performing the transformation on -ffcf_D_c2 
from atomic basis to the fragment basis (which are from 
the eignvectors of molecular fragment CI — U — C2 and 
are orthonormal), we take, 



lCB 



LOff 



-C2, m,n 



0, 



(rn = n), 
(m 7^ n), 



(22) 



with m, n are the index for the eigenvectors of the frag- 
ment CI — 



■ C2, and h^j^_ 



-C2 



are the the element of 



effective Hamiltonian matrix H°^. 

Within the fragment basis, from equation (|15p . we can 
obtain the self energy for CB correction in DFT as fol- 
lows. 
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FIG. 5: (Color online) The transmission function T{ui) in 
spectral space iv. The black solid curve is for the transmission 
function in the case without CB correction, and the red dash 
is with CB correction and Um = 1.0 eV, while the black dot- 
dash line is the Fermi level Ev = —8.5 eV. In the case without 
CB correction, we can see that there are seven levels close to 
the Fermi energy, the HOMO and the LUMO-l-1 are double- 
degenerate, respectively. When CB correction is introduced, 
for charging effect, the LUMO and the HOMO-2 will split, 
respectively. See text for details. 



rn.m.a \ m.a / 



where 



Urr 



fragment orbital m, and Sm.ri'o- 



- (1 - {rira.s)) ■ Um ' 

(23) 

is the Hubbard energy of the 
if m ^ n. The occu- 



can be ob- 



pation number of electrons {nm,s) = p. 
tained by the transformation on the density matrix from 
atomic basis Pfi^u.s to fragment basis Pm,n,s- Considering 
DCC along the idea similar to the case in LDA-|-Ur^i2£ 
we can get 6°, = /io-D-C2,m,m,a " U,n- The 

way to calculate Um within DFTB/gDFTB is shown in 
appendix IB] 

Finally the CB-correction self energy S^^^)'" will 



(CB)r 



be transformed back to the atomic basis E 

from the fragment basis Em,m,ff- By the help of 
Eqs. Iini), (0, (HZl) and HH]), introducing the self en- 
ergy I](")''/'^/< from leads, we can calculate the GFs as 
follows, 



"-^Cl-D-C2 ^ 1 ^Cl 



KS 



-C2 



and 



~C2 



Gr Y^(Q)</oa 



;(CB)r/ 



-C2- 



(24) 
(25) 



It should be noted that the eigenvectors of the frag- 
ment CI — D — C2 is updated in every cycle according to 
the updated KS Hamiltonian. 



B. calculation on a toy molecule 

The scheme above is realized within gBFTB^iSiiiS A 
toy molecule (C2S2) is taken for the calculation (see 
Fig. 1). The C2 is the CB part, while S is the linker. 
The bond length of C-C is 1.2 A, and the one of S-C is 
3.0 A. 

In the testing calculation of this chapter, for clarity 
within the level structure and the CB characteristics, the 
contribution from fragments CI and C2 is ignored, and 
Um = 1.0 eV is taken instead of the Um calculated by the 
method in appendix IB] (which are about 9.0 eV). Also for 
simplicity in the calculation, the fictitious golden leads 
are uscd)^ and the minimal basis is taken. 

The transmission function in spectral space is shown 
in Fig. 5. In the case without CB correction, we can see 
that there are seven levels close to the Fermi energy, while 
the HOMO and the LUMO+1 are double-degenerate, re- 
spectively. There is no level-split for charging, though 
the electronic occupation numbers of the seven levels 
above arc: 2.00, 1.89, 1.99, 1.99, 0.15, 0.00, 0.00, (it 
is clear that LUMO and HOMO-2 are not completely 
empty/occupied). Then, introducing CB correction, the 
LUMO and the HOMO-2 wiU split into em and + Um 
for the charging effect in the open shell, respectively. 

The electronic occupation number of LUMO and the 
current as a function of bias voltage Vbias are shown 
in Fig. 6. The charging-induced steps and the correct 
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be opened, which is the reason of the '1/3' in the ratio. 



IV. CONCLUSION 

In this paper, we have shown a new way to introduce 
CB correction to DFT calculation for the non-equilibrium 
case. The main elements of the approach are the follow- 
ing. 

1) With the help of the Ansatz in Ref J^, a self-energy 
is proposed for the (non-equilibrium) CB by the model 
Hamiltonian and EOM approach. From the comparison 
with the complete CB results, we can see that it can in- 
clude the characteristics of CB correctly. Further, the 
more important is that, from the view of DFT-based 
quantum transport calculation in non- equilibrium case, 
this self-energy is very convenient for programming more 
than the truncation in RefJ^ and ME approaches. 

2) Based on this self energy, a scheme with DCC is 
proposed to introduce CB correction to non- equilibrium 
DFT calculation. As is then reahzed within gOFTB^i^iiS 

By the new code above, the quantum-transport prop- 
erties of a toy molecule is calculated in CB regime. In the 
results of the electronic occupation number and the cur- 
rent as a function of bias voltage, the CB characteristics 
(the charging- induced steps and the ratio of 2/3 : 1/3 in 
the step heights) appear correctly. 



FIG. 6: (Color online) (a) the electronic occupation number 
T^LUMO of LUMO as a function of the bias voltage, (b) the 
current 7 as a function of the bias voltage. In them, the CB 
characteristics (the charging-induced steps and the ratio of 
2/3 : 1/3 in the heights of the steps) are very clear. 



ratio of 2/3 : 1/3 in the heights of the steps are very 
clear, which is consistent with the results in Fig. 2. 
The physics under it can be understood with the help 
of the charging-induced level-split, which is shown in 
the following. 1) When Hias = 0, the LUMO level 
is almost empty. 2) Then the positive bias voltage is 
added, and when the LUMO level elumo is coming into 
the Fermi windows (/ip — /Xp), there are two channels 
(elumct = ELUMoa) to be opened for current. As 
leads to the '2/3' in the ratio. 3) After that, the elec- 
tronic occupation number of LUMO comes to be 0.66 
(see Fig. 6(a) and Fig. 2(a)), and for the charging effect, 
the degenerate levels (clumcTj ^LUMO.i) will split into 
two (eLUMO,T and eLUMO.i = eLUMO,T + t^LUMo)- 4) At 
the time that the level eLUMO,i(= elumoj + f^LUMo) en- 
ters the Fermi windows, there will be only one channel to 
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APPENDIX A: WAY TO CALCULATE e'^'^X 

For simplicity, here, only a double-level (m = 1 with 
spin-up and spin-down) case is taken to describe the 
quantum dot. Then the index m is ignored in this part, 
and from Eq. ([2]), the Hamiltonian of the dot can be re- 
written as, 

Hu=Y.eyU^ + lY.Un,n^. (Al) 

(7 (7 
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By Eq. (ITT|), with the help of the Ansatz in RcfiiS, we 

can directly write out G^^^^ as follows, 



LO — e'^ — U — h] 



2^i/H<5(c.-e|^)(l-(,^^)) 
+2TTiJ{uj)6{u-el-U){n^) 
2i7?7(a;)(l-(n^)) , 2ii^J{uj){ns) 



(A2) 



where /(w) is sojne fcinrf 0/ pseudo-Fermi function in non- 
equilibrium case, and 77 0+ 

Assuming that there exists the relation E^*^^'^ ~ 
i/(ijj)r('-'^' (ti>), and G'^'^^ can be re-written in the form. 



QiUX ^ Gi^)'-s(CB)<G^c/)a^ (A3) 
Then with the help of Eq. PT]) , we will obtain 



(A4) 



Comparing Eqs. (jA2p and (|A4p . we can get 

T^'^^\uj) = 2i^D{uj), (A5) 

with 

1 " {n-a) {ns)_ 

i^(co) = ^ ^ 7 ^ ^ " . (A6) 



(1 - {n,)r 



Therefore, it is clear that si™'^ = if{uj)T'^'^^^Lu) 0, 
since 77 ^ 0^ and D{uj) never comes to be infinite. 



APPENDIX B: WAY TO CALCULATE 
HUBBARD ENERGY U BY DFTB/GDFTB 

In DFTB/gDFTB^i^iiaiiii22i33 ^^le tight-binding 
Hamiltonian is written as follows, 



where fx G a and ly G f3. i is the index for the eigen- 
state of molecule (or molecular fragment), fx^v are the 
index of atomic basis, while a, (3, ^ indicates atoms, qa is 
the charge of atom a, and rii is the electron occupation 
number of eigenstate i. "/a,f3 is a function of Ua, Up and 
|R„ - R/3|,-^'^ with Ua {Up) is the Hubbard U of atom 
a (atom /?), and |Rq — R^| is the distance between the 
two atoms, c^.i is the project of the eigenvector i on the 
atomic basis /i, while S'^.j, are the elements of overlap 
matrix in atomic basis. 

With the help of linear combination of atomic or- 
bitals(LCAO) ansatz,— we can transform the Hamilto- 
nian matrix (jBip from atomic basis to molecular basis, 
and get the eigenvalue of the eigenstate i as follows, 



(B4) 



Then, according to the defination of Hubbard U )35,36,37 
we will obtain, 



dtij dtij drij duj 



sec 



E 



dc 



,dc„ 



dn 



(B5) 



Ignoring the contribution from dc/ dn and dH'^ /dn, with 
the condition dni/dnj = Sij, we can get 



dn, '^'^ 



N 



N 



with 



TT ttO I rrSCC 

N 



(Bl) 



H^T = j5^,.E(7a,?+7/3.«)Ag« (B2) 



^ occ 



with fx G a, V G f3. 

The calculation of Hubbard [/ by the approxima- 
tions (jB6P is performed in four examples: 1) guanine- 
cytosine base pair (GC), 2) adenine-thymine base pair 
(AT), 3) benzene (CeHg), 4) double carbon (C2) in the 
case that the bond length is 1.2 A. The results are shown 
in table I. The error from the approximations is less than 
12%, which is acceptable within DFTB/gDFTB. 
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TABLE I: The calculation of Hubbard U by approxima- 
tion (|B6)) 





GC 


AT 




C2 


C/homo (eV) ^ 


6.170 


6.235 


7.157 


9.040 


deaoMO (gV) 


5.556 


5.645 


6.451 


9.022 


Ai (eV) 


0.614 


0.590 


0.706 


0.018 


A2 (%) ^ 


11.05 


10.45 


10.94 


0.20 



^ results of HOMO by approximation IIB6I 1. 
Ai = J7homo - 9<:homo/9"homo- 

" A2 = Ai/(SeH0M0/9"H0M0)- 
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